#################################################
# Code to replicate monadic analysis for:       #
# "Civil War Mediation and Integration into     #
# Global Value Chains"                          #
#                                               #
# J. Tyson Chatagnier                           #
#                                               #
# To be published in International Interactions #
#                                               #
# File created 16 September 2018                #
#################################################

rm(list = ls())

library(readstata13)
library(ggplot2)
library(texreg)
library(countrycode)

#######################
# Function to compute #
# substantive effects #
####################### 

AvgPreds <- function(
                     mod,
                     xmat       = NULL,
                     varcov     = NULL,
                     var.change,
                     vars       = NULL,
                     var.hold   = NULL,
                     change.lo  = NULL,
                     change.hi  = NULL,
                     change.vec = NULL,
                     n.out      = 1000
                     ) {
    if (is.null(xmat) & !is.numeric(mod$x)) {
        cat("Please provide an x matrix.\n")
    } else if (is.null(xmat) & is.numeric(mod$x)) {
        xmat <- mod$x
    }
    k <- ncol(xmat)
    b <- mod$coefficients
    if (is.null(change.vec)) {
        change.lo <- ifelse(
                            is.null(change.lo),
                            min(xmat[, var.change]),
                            change.lo
                            )
        change.hi <- ifelse(
                            is.null(change.hi),
                            max(xmat[, var.change]),
                            change.hi
                            )
        change.vec <- seq(
                          from       = change.lo,
                          to         = change.hi,
                          length.out = n.out
                          )
    }
    xc <- xmat
    if (!is.null(var.hold)) {
      xc[, var.hold[1]] <- var.hold[2]
    }
    pred.mat <- matrix(
                       nrow = length(change.vec),
                       ncol = 3
                       )
    for (i in 1 : length(change.vec)) {
      xc[, var.change] <- change.vec[i]
      if (!is.null(vars)) {
        for (j in 1 : nrow(vars)) {
          xc[, vars[j, 1]] <- vars[j, 2]
        }
                        # vars is an optional variable that takes
                        # a matrix as its argument, with the first
                        # column being variable locations to be held
                        # constant and the second being the values
                        # at which to hold those variables.
      }
      pp <- plogis(xc %*% b)    
      deriv <- dlogis(xc %*% b)
                      # dlogis() is the derivative of plogis()
      deriv.mat <- sweep(
                         xc,
                         1,
                         deriv,
                         '*'
                         )
                        # Multiply each row of xc by each value of deriv.
      if (is.null(varcov)) {
        vhat <- vcov(mod)
      } else {
        vhat <- varcov
      }
      se.pp <- sqrt(rowSums((deriv.mat %*% vhat) * deriv.mat))
                        # Just need the standard errors, so this is
                        # quicker than matrix multiplication for the
                        # entire covariance matrix, and it avoids 
                        # problems with very large data sets.
      ci.lo <- pp - (1.96 * se.pp)
      ci.lo <- ifelse(
                      ci.lo < 0,
                      0,
                      ifelse(
                             ci.lo > 1,
                             1,
                             ci.lo
                             )
                      )
      ci.hi <- pp + (1.96 * se.pp)
      ci.hi <- ifelse(
                      ci.hi < 0,
                      0,
                      ifelse(
                             ci.hi > 1,
                             1,
                             ci.hi
                             )
                      )
      pred.mat[i, 1] <- mean(pp)
      pred.mat[i, 2] <- mean(ci.lo)
      pred.mat[i, 3] <- mean(ci.hi)
    }
  return(
         list(
              pred   = pred.mat[, 1],
              lo     = pred.mat[, 2],
              hi     = pred.mat[, 3],
              change = change.vec
              )
         )
}

#########################
# Read in the data sets #
#########################

 # setwd("")
                        # Fill in the area between the quotes and uncomment
                        # to specify a working directory.

data.set <- read.dta13("gvc_mediation_replication_data.dta")

####################
# Run Logit models #
####################

mod.val <- glm(
	             majpow_med ~
	                          constpc +
	                          ln.tot +
	                          colonial +
	                          polity2 +
	                          cinc + 
	                          africa +
	                          asia +
                            num.fail +
                            num.ongo +
	                          ln.exp + 
	                          t +
	                          t2 +
	                          t3,
	             data   = data.set,
	             family = binomial(link = "logit"),
	             x      = TRUE
	             )

mod.glob <- glm(
	              majpow_med ~
	                           constpc +
	                           ln.tot +
	                           colonial +
	                           polity2 +
	                           cinc + 
	                           africa +
	                           asia +
                             num.fail +
                             num.ongo +
	                           expshare.glob + 
	                           t +
	                           t2 +
	                           t3,
	              data   = data.set,
	              family = binomial(link = "logit"),
	              x      = TRUE
	              )

val.imp <- glm(
	             majpow_med ~
	                          constpc +
	                          ln.tot +
	                          colonial +
	                          polity2 +
	                          cinc + 
	                          africa +
	                          asia +
                            num.fail +
                            num.ongo +
	                          ln.exp + 
	                          ln.imp +
	                          t +
	                          t2 +
	                          t3,
	             data   = data.set,
	             family = binomial(link = "logit"),
	             x      = TRUE
	             )

glob.imp <- glm(
	              majpow_med ~
	                           constpc +
	                           ln.tot +
	                           colonial +
	                           polity2 +
	                           cinc + 
	                           africa +
	                           asia +
                             num.fail +
                             num.ongo +
	                           expshare.glob + 
	                           impshare.glob +
	                           t +
	                           t2 +
	                           t3,
	              data   = data.set,
	              family = binomial(link = "logit"),
	              x      = TRUE
	              )

############
# Table 1  #
############

screenreg(
          list(
               mod.val,
               mod.glob,
               val.imp,
               glob.imp
               ),
          stars              = c(
                                 0.01,
                                 0.05,
                                 0.1
                                 ) ,
          custom.coef.names  = c(
                                 "Intercept",
                                 "GDP per capita",
                                 "Total trade (logged)",
                                 "Colonial ties",
                                 "Polity score",
                                 "Capabilities",
                                 "Africa",
                                 "Asia",
                                 "Failed mediation",
                                 "Ongoing mediation",
                                 "Intermediate exports (logged)",
                                 "Time",
                                 "Time$^{2}$",
                                 "Time$^{3}$",
                                 "Share of intermediate exports",
                                 "Intermediate imports (logged)",
                                 "Share of intermediate imports"
                                 ),
          custom.model.names = c(
                                 "Logged Value",
                                 "Global Share",
                                 "Logged Value",
                                 "Global Share"
                                 )
          )

#######################
# Substantive effects #
#######################

n.vals <- 1000

imp.lo <- mean(data.set$impshare.glob) - (sd(data.set$impshare.glob) / 2)
imp.hi <- mean(data.set$impshare.glob) + (sd(data.set$impshare.glob) / 2)

exp.lo <- mean(data.set$expshare.glob) - (sd(data.set$expshare.glob) / 2)
exp.hi <- mean(data.set$expshare.glob) + (sd(data.set$expshare.glob) / 2)
                        # We wants shifts of 1/2 of a standard deviation
                        # above and below the mean for our variables of
                        # interest. Create those here.

###################################
# Get average predicted values    #
# as we vary intermediate         #
# exports about its mean, holding #
# intermediate imports constant   #
# at three different values       #
###################################

exp.lowimp <- AvgPreds(
                       mod        = glob.imp,
                       var.change = 11,
                       var.hold   = c(
                                      12,
                                      imp.lo
                                      ),
                       n.out      = n.vals
                       )

exp.medimp <- AvgPreds(
                       mod        = glob.imp,
                       var.change = 11,
                       var.hold   = c(
                                      12,
                                      mean(data.set$impshare.glob)
                                      ),
                       n.out      = n.vals
                       )

exp.highimp <- AvgPreds(
                        mod        = glob.imp,
                        var.change = 11,
                        var.hold   = c(
                                       12,
                                       imp.hi
                                       ),
                        n.out      = n.vals
                        )

###################################
# Get average predicted values    #
# as we vary intermediate         #
# imports about its mean, holding #
# intermediate exports constant   #
# at three different values       #
###################################

imp.lowexp <- AvgPreds(
                       mod        = glob.imp,
                       var.change = 12,
                       var.hold   = c(
                                      11,
                                      exp.lo
                                      ),
                       n.out      = n.vals
                       )

imp.medexp <- AvgPreds(
                       mod        = glob.imp,
                       var.change = 12,
                       var.hold   = c(
                                      11,
                                      mean(data.set$expshare.glob)
                                      ),
                       n.out      = n.vals
                       )

imp.highexp <- AvgPreds(
                        mod        = glob.imp,
                        var.change = 12,
                        var.hold   = c(
                                       11,
                                       exp.hi
                                       ),
                        n.out      = n.vals
                        )

###########################
# Create data frames with #
# results from predicted  #
# probabilities in each   #
# case                    #
###########################

exp.df <- data.frame(
	                   exp = c(
	                   	       exp.lowimp$change,
	                   	       exp.medimp$change,
                             exp.highimp$change
	                   	       ),
	                   p   = c(
	                   	       exp.lowimp$pred,
                             exp.medimp$pred,
	                   	       exp.highimp$pred
	                   	       ),
	                   lo  = c(
	                   	       exp.lowimp$lo,
                             exp.medimp$lo,
	                   	       exp.highimp$lo
	                   	       ),
	                   hi  = c(
	                   	       exp.lowimp$hi,
                             exp.medimp$hi,
	                   	       exp.highimp$hi
	                   	       ),
	                   mod = rep(
	                   	         c(
	                   	         	 "Low",
	                   	         	 "Mean",
                                 "High"
	                   	         	 ),
	                   	         each = n.vals
	                   	         )
	                   )

exp.df$mod <- factor(
                     exp.df$mod, 
                     levels = c(
                                "Low",
                                "Mean",
                                "High"
                                )
                     )

imp.df <- data.frame(
                     imp = c(
                             imp.lowexp$change,
                             imp.medexp$change,
                             imp.highexp$change
                             ),
                     p   = c(
                             imp.lowexp$pred,
                             imp.medexp$pred,
                             imp.highexp$pred
                             ),
                     lo  = c(
                             imp.lowexp$lo,
                             imp.medexp$lo,
                             imp.highexp$lo
                             ),
                     hi  = c(
                             imp.lowexp$hi,
                             imp.medexp$hi,
                             imp.highexp$hi
                             ),
                     mod = rep(
                               c(
                                 "Low",
                                 "Mean",
                                 "High"
                                 ),
                               each = n.vals
                               )
                     )

imp.df$mod <- factor(
                     imp.df$mod, 
                     levels = c(
                                "Low",
                                "Mean",
                                "High"
                                )
                     )

##########################
# Left panel of Figure 1 #
##########################

ggplot(
       data = exp.df,
       aes(
           x        = exp,
           y        = p,
           colour   = mod,
           fill     = mod,
           linetype = mod
           )
       ) +
  geom_ribbon(
              aes(
                  ymin = lo,
                  ymax = hi
                  ),
              alpha  = 0.3,
              colour = NA
              ) +
  geom_line(size   = 1.10) +
  geom_rug(
           data        = data.set, 
           aes(
               x = expshare.glob
               ), 
           sides       = "b", 
           inherit.aes = F, 
           colour      = "navyblue", 
           size        = 0.5
           ) +
  labs(
       fill     = "Import Level", 
       colour   = "Import Level", 
       linetype = "Import Level"
       ) +
  xlab("Percent Share of Global Intermediate Exports") +
  ylab("Predicted Probability of Mediation") +
  theme_bw()

###########################
# Right panel of Figure 1 #
###########################

ggplot(
       data = imp.df,
       aes(
           x        = imp,
           y        = p,
           colour   = mod,
           fill     = mod,
           linetype = mod
           )
       ) +
  geom_ribbon(
              aes(
                  ymin = lo,
                  ymax = hi
                  ),
              alpha  = 0.3,
              colour = NA
              ) +
  geom_line(size   = 1.10) +
  geom_rug(
           data        = data.set, 
           aes(
               x = impshare.glob
               ), 
           sides       = "b", 
           inherit.aes = F, 
           colour      = "navyblue", 
           size        = 0.5
           ) +
  labs(
       fill     = "Export Level", 
       colour   = "Export Level", 
       linetype = "Export Level"
       ) +
  xlab("Percent Share of Global Intermediate Imports") +
  ylab("Predicted Probability of Mediation") +
  theme_bw()

##################################
#                                #
#            APPENDIX            #
#                                #
##################################

######################
# Summary statistics #
######################

varnames <- c(
              "Mediation onset (all)",
              "Mediation onset (major powers)",
              "Intermediate exports (logged)",
              "Intermediate imports (logged)",
              "Intermediate exports (share)",
              "Intermediate imports (share)",
              "GDP per capita",
              "Total trade (logged)",
              "Colonial ties",
              "Polity score",
              "Capabilities",
              "Africa",
              "Asia",
              "Failed mediations",
              "Ongoing mediations",
              "Year",
              "Years since conflict onset"
              )

vars <- cbind(
              data.set$Med_yes.no,
              data.set$majpow_med,
              data.set$ln.exp,
              data.set$ln.imp,
              data.set$expshare.glob,
              data.set$impshare.glob,
              data.set$constpc,
              data.set$ln.tot,
              data.set$colonial,
              data.set$polity2,
              data.set$cinc,
              data.set$africa,
              data.set$asia,
              data.set$num.fail,
              data.set$num.ongo,
              data.set$year,
              data.set$t
              )

n.vars <- apply(
                vars,
                2,
                function (x) {
                  length(na.omit(x))
                }
                )

mean.vars <- apply(
                   vars,
                   2,
                   mean,
                   na.rm = TRUE
                   )

mean.vars[16] <- trunc(mean.vars[16])
                        # Truncate the year variable so 
                        # that it's a whole number

sd.vars <- apply(
                 vars,
                 2,
                 sd,
                 na.rm = TRUE
                 )

min.vars <- apply(
                  vars,
                  2,
                  min,
                  na.rm = TRUE
                  )

max.vars <- apply(
                  vars,
                  2,
                  max,
                  na.rm = TRUE
                  )


summary.df <- data.frame(
                         Variable     = varnames,
                         Observations = n.vars,
                         Mean         = mean.vars,
                         SD           = sd.vars,
                         Min          = min.vars,
                         Max          = max.vars
                         )

summary.df[, 3 : 6] <- round(
                             summary.df[, 3 : 6], 
                             3
                             )

############
# Table A1 #
############

print(summary.df)

#################################
# Analysis using OECD countries #
#################################

#####################
# Read in OECD data #
#####################

data.oecd <- read.dta13("gvc_mediation_replication_data_oecd.dta")

##############################
# Run analyses using OECD as #
# dependent variable         #
##############################

mod.val.oecd <- glm(
                    oecd_med ~
                                 constpc +
                                 ln.tot +
                                 colonial +
                                 polity2 +
                                 cinc + 
                                 africa +
                                 asia +
                                 num.fail +
                                 num.ongo +
                                 ln.exp + 
                                 # ln.imp +
                                 t +
                                 t2 +
                                 t3,
                    data   = data.oecd,
                    family = binomial(link = "logit"),
                    x      = TRUE
                    )

mod.glob.oecd <- glm(
                     oecd_med ~
                                  constpc +
                                  ln.tot +
                                  colonial +
                                  polity2 +
                                  cinc + 
                                  africa +
                                  asia +
                                  num.fail +
                                  num.ongo +
                                  expshare.glob + 
                                  # impshare.glob +
                                  t +
                                  t2 +
                                  t3,
                     data   = data.oecd,
                     family = binomial(link = "logit"),
                     x      = TRUE
                     )

val.imp.oecd <- glm(
                    oecd_med ~
                                 constpc +
                                 ln.tot +
                                 colonial +
                                 polity2 +
                                 cinc + 
                                 africa +
                                 asia +
                                 num.fail +
                                 num.ongo +
                                 ln.exp + 
                                 ln.imp +
                                 t +
                                 t2 +
                                 t3,
                    data   = data.oecd,
                    family = binomial(link = "logit"),
                    x      = TRUE
                    )

glob.imp.oecd <- glm(
                     oecd_med ~
                                  constpc +
                                  ln.tot +
                                  colonial +
                                  polity2 +
                                  cinc + 
                                  africa +
                                  asia +
                                  num.fail +
                                  num.ongo +
                                  expshare.glob + 
                                  impshare.glob +
                                  t +
                                  t2 +
                                  t3,
                     data   = data.oecd,
                     family = binomial(link = "logit"),
                     x      = TRUE
                     )

############
# Table A2 #
############

screenreg(
          list(
               mod.val.oecd,
               mod.glob.oecd,
               val.imp.oecd,
               glob.imp.oecd
               ),
          stars              = c(
                                 0.01,
                                 0.05,
                                 0.1
                                 ) ,
          custom.coef.names  = c(
                                 "Intercept",
                                 "GDP per capita",
                                 "Total trade (logged)",
                                 "Colonial ties",
                                 "Polity score",
                                 "Capabilities",
                                 "Africa",
                                 "Asia",
                                 "Failed mediations",
                                 "Ongoing mediations",
                                 "Intermediate exports (logged)",
                                 "Time",
                                 "Time$^{2}$",
                                 "Time$^{3}$",
                                 "Share of intermediate exports",
                                 "Intermediate imports (logged)",
                                 "Share of intermediate imports"
                                 ),
          custom.model.names = c(
                                 "Logged Value",
                                 "Global Share",
                                 "Logged Value",
                                 "Global Share"
                                 )
          )

############################
# Excluding potentially    #
# influential observations #
############################

################################
# Create plot of import values #
################################

data.set$abbrev <- countrycode(
                               data.set$ccode,
                               "cown",
                               "cowc"
                               )
                        # Get abbreviations for plot

avg.imp <- cbind(
                 by(
                    data.set$ln.imp,
                    list(data.set$abbrev),
                    mean
                    )
                 )
                        # Find average intermediate
                        # imports for each country

num.obs <- cbind(
                 by(
                    data.set$ln.imp,
                    list(data.set$abbrev),
                    length
                    )
                 )
                        # Find how many times each country
                        # appears in the data

impval.df <- data.frame(
                        imports = avg.imp,
                        obs     = num.obs,
                        abbrev  = rownames(num.obs)
                        )

#############
# Figure A1 #
#############

ggplot(
       impval.df, 
       aes_string(
                  x = "imports", 
                  y = "obs"
                  )
                  ) + 
         geom_text(
                   size  = 3, 
                   aes(label = abbrev)
                   ) +
         xlab("Average intermediate imports (logged)") +
         ylab("Number of dispute-years") +
         theme_bw()

############################
# Omit India from analysis #
############################

mod.val.noind <- glm(
                     majpow_med ~
                                  constpc +
                                  ln.tot +
                                  colonial +
                                  polity2 +
                                  cinc + 
                                  africa +
                                  asia +
                                  num.fail +
                                  num.ongo +
                                  ln.exp + 
                                  # ln.imp +
                                  t +
                                  t2 +
                                  t3,
                     data   = data.set,
                     subset = abbrev != "IND",
                     family = binomial(link = "logit"),
                     x      = TRUE
                     )

mod.glob.noind <- glm(
                      majpow_med ~
                                   constpc +
                                   ln.tot +
                                   colonial +
                                   polity2 +
                                   cinc + 
                                   africa +
                                   asia +
                                   num.fail +
                                   num.ongo +
                                   expshare.glob + 
                                   # impshare.glob +
                                   t +
                                   t2 +
                                   t3,
                      data   = data.set,
                      subset = abbrev != "IND",
                      family = binomial(link = "logit"),
                      x      = TRUE
                      )

val.imp.noind <- glm(
                     majpow_med ~
                                  constpc +
                                  ln.tot +
                                  colonial +
                                  polity2 +
                                  cinc + 
                                  africa +
                                  asia +
                                  num.fail +
                                  num.ongo +
                                  ln.exp + 
                                  ln.imp +
                                  t +
                                  t2 +
                                  t3,
                     data   = data.set,
                     subset = abbrev != "IND",
                     family = binomial(link = "logit"),
                     x      = TRUE
                     )

glob.imp.noind <- glm(
                      majpow_med ~
                                   constpc +
                                   ln.tot +
                                   colonial +
                                   polity2 +
                                   cinc + 
                                   africa +
                                   asia +
                                   num.fail +
                                   num.ongo +
                                   expshare.glob + 
                                   impshare.glob +
                                   t +
                                   t2 +
                                   t3,
                      data   = data.set,
                      subset = abbrev != "IND",
                      family = binomial(link = "logit"),
                      x      = TRUE
                      )

############
# Table A3 #
############
      
screenreg(
          list(
               mod.val.noind,
               mod.glob.noind,
               val.imp.noind,
               glob.imp.noind
               ),
          stars              = c(
                                 0.01,
                                 0.05,
                                 0.1
                                 ) ,
          custom.coef.names  = c(
                                 "Intercept",
                                 "GDP per capita",
                                 "Total trade (logged)",
                                 "Colonial ties",
                                 "Polity score",
                                 "Capabilities",
                                 "Africa",
                                 "Asia",
                                 "Failed mediation",
                                 "Ongoing mediation",
                                 "Intermediate exports (logged)",
                                 "Time",
                                 "Time$^{2}$",
                                 "Time$^{3}$",
                                 "Share of intermediate exports",
                                 "Intermediate imports (logged)",
                                 "Share of intermediate imports"
                                 ),
          custom.model.names = c(
                                 "Logged Value",
                                 "Global Share",
                                 "Logged Value",
                                 "Global Share"
                                 )
          )

#######################################
# Omit major intermediate importers   #
# from the analysis                   #
#######################################

mod.val.implo <- glm(
                     majpow_med ~
                                  constpc +
                                  ln.tot +
                                  colonial +
                                  polity2 +
                                  cinc + 
                                  africa +
                                  asia +
                                  num.fail +
                                  ln.exp + 
                                  t +
                                  t2 +
                                  t3,
                     data   = data.set,
                     subset = ln.imp < 15,
                     family = binomial(link = "logit"),
                     x      = TRUE
                     )

mod.glob.implo <- glm(
                      majpow_med ~
                                   constpc +
                                   ln.tot +
                                   colonial +
                                   polity2 +
                                   cinc + 
                                   africa +
                                   asia +
                                   num.fail +
                                   expshare.glob + 
                                   t +
                                   t2 +
                                   t3,
                      data   = data.set,
                      subset = ln.imp < 15,
                      family = binomial(link = "logit"),
                      x      = TRUE
                      )

val.imp.implo <- glm(
                     majpow_med ~
                                  constpc +
                                  ln.tot +
                                  colonial +
                                  polity2 +
                                  cinc + 
                                  africa +
                                  asia +
                                  num.fail +
                                  ln.exp + 
                                  ln.imp +
                                  t +
                                  t2 +
                                  t3,
                     data   = data.set,
                     subset = ln.imp < 15,
                     family = binomial(link = "logit"),
                     x      = TRUE
                     )

glob.imp.implo <- glm(
                      majpow_med ~
                                   constpc +
                                   ln.tot +
                                   colonial +
                                   polity2 +
                                   cinc + 
                                   africa +
                                   asia +
                                   num.fail +
                                   expshare.glob + 
                                   impshare.glob +
                                   t +
                                   t2 +
                                   t3,
                      data   = data.set,
                      subset = ln.imp < 15,
                      family = binomial(link = "logit"),
                      x      = TRUE
                      )
      
############
# Table A4 #
############

screenreg(
          list(
               mod.val.implo,
               mod.glob.implo,
               val.imp.implo,
               glob.imp.implo
               ),
          stars              = c(
                                 0.01,
                                 0.05,
                                 0.1
                                 ) ,
          custom.coef.names  = c(
                                 "Intercept",
                                 "GDP per capita",
                                 "Total trade (logged)",
                                 "Colonial ties",
                                 "Polity score",
                                 "Capabilities",
                                 "Africa",
                                 "Asia",
                                 "Failed mediation",
                                 "Intermediate exports (logged)",
                                 "Time",
                                 "Time$^{2}$",
                                 "Time$^{3}$",
                                 "Share of intermediate exports",
                                 "Intermediate imports (logged)",
                                 "Share of intermediate imports"
                                 ),
          custom.model.names = c(
                                 "Logged Value",
                                 "Global Share",
                                 "Logged Value",
                                 "Global Share"
                                 )
          )

##########################
# Omit major powers from #
# the analysis           #
##########################

mp <- c(
        2,
        200,
        220,
        255,
        365,
        710,
        740
        )
                        # Set of major power countries
                        # for the time period in this
                        # analysis

mod.val.nomp <- glm(
                    majpow_med ~
                                 constpc +
                                 ln.tot +
                                 colonial +
                                 polity2 +
                                 cinc + 
                                 africa +
                                 asia +
                                 num.fail +
                                 num.ongo +
                                 ln.exp + 
                                 t +
                                 t2 +
                                 t3,
                    data   = data.set,
                    subset = !(ccode %in% mp),
                    family = binomial(link = "logit"),
                    x      = TRUE
                    )

mod.glob.nomp <- glm(
                     majpow_med ~
                                  constpc +
                                  ln.tot +
                                  colonial +
                                  polity2 +
                                  cinc + 
                                  africa +
                                  asia +
                                  num.fail +
                                  num.ongo +
                                  expshare.glob + 
                                  t +
                                  t2 +
                                  t3,
                     data   = data.set,
                     subset = !(ccode %in% mp),
                     family = binomial(link = "logit"),
                     x      = TRUE
                     )

val.imp.nomp <- glm(
                    majpow_med ~
                                 constpc +
                                 ln.tot +
                                 colonial +
                                 polity2 +
                                 cinc + 
                                 africa +
                                 asia +
                                 num.fail +
                                 num.ongo +
                                 ln.exp + 
                                 ln.imp +
                                 t +
                                 t2 +
                                 t3,
                    data   = data.set,
                    subset = !(ccode %in% mp),
                    family = binomial(link = "logit"),
                    x      = TRUE
                    )

glob.imp.nomp <- glm(
                     majpow_med ~
                                  constpc +
                                  ln.tot +
                                  colonial +
                                  polity2 +
                                  cinc + 
                                  africa +
                                  asia +
                                  num.fail +
                                  num.ongo +
                                  expshare.glob + 
                                  impshare.glob +
                                  t +
                                  t2 +
                                  t3,
                     data   = data.set,
                     subset = !(ccode %in% mp),
                     family = binomial(link = "logit"),
                     x      = TRUE
                     )
      
############
# Table A5 #
############

screenreg(
          list(
               mod.val.nomp,
               mod.glob.nomp,
               val.imp.nomp,
               glob.imp.nomp
               ),
          stars              = c(
                                 0.01,
                                 0.05,
                                 0.1
                                 ) ,
          custom.coef.names  = c(
                                 "Intercept",
                                 "GDP per capita",
                                 "Total trade (logged)",
                                 "Colonial ties",
                                 "Polity score",
                                 "Capabilities",
                                 "Africa",
                                 "Asia",
                                 "Failed mediation",
                                 "Ongoing mediation",
                                 "Intermediate exports (logged)",
                                 "Time",
                                 "Time$^{2}$",
                                 "Time$^{3}$",
                                 "Share of intermediate exports",
                                 "Intermediate imports (logged)",
                                 "Share of intermediate imports"
                                 ),
          custom.model.names = c(
                                 "Logged Value",
                                 "Global Share",
                                 "Logged Value",
                                 "Global Share"
                                 )
          )

############################
# Mediation by any country #
############################

mod.val.any <- glm(
                   Med_yes.no ~
                                constpc +
                                ln.tot +
                                colonial +
                                polity2 +
                                cinc + 
                                africa +
                                asia +
                                ln.exp + 
                                t +
                                t2 +
                                t3,
                   data   = data.set,
                   family = binomial(link = "logit"),
                   x      = TRUE
                   )

mod.glob.any <- glm(
                    Med_yes.no ~
                                 constpc +
                                 ln.tot +
                                 colonial +
                                 polity2 +
                                 cinc + 
                                 africa +
                                 asia +
                                 expshare.glob + 
                                 t +
                                 t2 +
                                 t3,
                    data   = data.set,
                    family = binomial(link = "logit"),
                    x      = TRUE
                    )

val.imp.any <- glm(
                   Med_yes.no ~
                                constpc +
                                ln.tot +
                                colonial +
                                polity2 +
                                cinc + 
                                africa +
                                asia +
                                ln.exp + 
                                ln.imp +
                                t +
                                t2 +
                                t3,
                   data   = data.set,
                   family = binomial(link = "logit"),
                   x      = TRUE
                   )

glob.imp.any <- glm(
                    Med_yes.no ~
                                  constpc +
                                  ln.tot +
                                  colonial +
                                  polity2 +
                                  cinc + 
                                  africa +
                                  asia +
                                  expshare.glob + 
                                  impshare.glob +
                                  t +
                                  t2 +
                                  t3,
                     data   = data.set,
                     family = binomial(link = "logit"),
                     x      = TRUE
                     )

############
# Table A6 #
############

screenreg(
          list(
               mod.val.any,
               mod.glob.any,
               val.imp.any,
               glob.imp.any
               ),
          stars              = c(
                                 0.01,
                                 0.05,
                                 0.1
                                 ) ,
          custom.coef.names  = c(
                                 "Intercept",
                                 "GDP per capita",
                                 "Total trade (logged)",
                                 "Colonial ties",
                                 "Polity score",
                                 "Capabilities",
                                 "Africa",
                                 "Asia",
                                 "Intermediate exports (logged)",
                                 "Time",
                                 "Time$^{2}$",
                                 "Time$^{3}$",
                                 "Share of intermediate exports",
                                 "Intermediate imports (logged)",
                                 "Share of intermediate imports"
                                 ),
          custom.model.names = c(
                                 "Logged Value",
                                 "Global Share",
                                 "Logged Value",
                                 "Global Share"
                                 )
          )

###########################
# Get substantive effects #
###########################

exp.medimp.any <- AvgPreds(
                           mod        = glob.imp.any,
                           var.change = 9,
                           var.hold   = c(
                                          10,
                                          mean(data.set$impshare.glob)
                                          ),
                           n.out      = n.vals
                           )

imp.medexp.any <- AvgPreds(
                           mod        = glob.imp.any,
                           var.change = 10,
                           var.hold   = c(
                                          9,
                                          mean(data.set$expshare.glob)
                                          ),
                           n.out      = n.vals
                           )

###############################
# Create relevant data frames #
###############################

exp.df.any <- data.frame(
                         exp = exp.medimp.any$change,
                         p   = exp.medimp.any$pred,
                         lo  = exp.medimp.any$lo,
                         hi  = exp.medimp.any$hi
                         )


imp.df.any <- data.frame(
                         imp = imp.medexp.any$change,
                         p   = imp.medexp.any$pred,
                         lo  = imp.medexp.any$lo,
                         hi  = imp.medexp.any$hi
                         )

###########################
# Left panel of Figure A2 #
###########################

ggplot(
       data = exp.df.any,
       aes(
           x        = exp,
           y        = p
           )
       ) +
  geom_ribbon(
              aes(
                  ymin = lo,
                  ymax = hi
                  ),
              alpha  = 0.3,
              colour = NA,
              fill   = "navyblue"
              ) +
  geom_line(
            size   = 1.10,
            colour = "navyblue"
            ) +
  xlab("Percent Share of Global Intermediate Exports") +
  ylab("Predicted Probability of Mediation") +
  theme_bw()

############################
# Right panel of Figure A2 #
############################

ggplot(
       data = imp.df.any,
       aes(
           x        = imp,
           y        = p
           )
       ) +
  geom_ribbon(
              aes(
                  ymin = lo,
                  ymax = hi
                  ),
              alpha  = 0.3,
              colour = NA,
              fill   = "navyblue"
              ) +
  geom_line(
            size   = 1.10,
            colour = "navyblue"
            ) +
  labs(
       fill     = "Export Level", 
       colour   = "Export Level", 
       linetype = "Export Level"
       ) +
  xlab("Percent Share of Global Intermediate Imports") +
  ylab("Predicted Probability of Mediation") +
  theme_bw()

###########################################
# Models without candidate-level controls #
###########################################

##############################
# Omit regime type and GDPpc #
##############################

mod.val.exclude <- glm(
                       majpow_med ~
                                    ln.tot +
                                    colonial +
                                    cinc + 
                                    africa +
                                    asia +
                                    num.fail +
                                    num.ongo +
                                    ln.exp + 
                                    t +
                                    t2 +
                                    t3,
                       data   = data.set,
                       family = binomial(link = "logit"),
                       x      = TRUE
                       )

mod.glob.exclude <- glm(
                        majpow_med ~
                                     ln.tot +
                                     colonial +
                                     cinc + 
                                     africa +
                                     asia +
                                     num.fail +
                                     num.ongo +
                                     expshare.glob + 
                                     t +
                                     t2 +
                                     t3,
                        data   = data.set,
                        family = binomial(link = "logit"),
                        x      = TRUE
                        )

val.imp.exclude <- glm(
                       majpow_med ~
                                    ln.tot +
                                    colonial +
                                    cinc + 
                                    africa +
                                    asia +
                                    num.fail +
                                    num.ongo +
                                    ln.exp + 
                                    ln.imp +
                                    t +
                                    t2 +
                                    t3,
                       data   = data.set,
                       family = binomial(link = "logit"),
                       x      = TRUE
                       )

glob.imp.exclude <- glm(
                        majpow_med ~
                                     ln.tot +
                                     colonial +
                                     cinc + 
                                     africa +
                                     asia +
                                     num.fail +
                                     num.ongo +
                                     expshare.glob + 
                                     impshare.glob +
                                     t +
                                     t2 +
                                     t3,
                        data   = data.set,
                        family = binomial(link = "logit"),
                        x      = TRUE
                        )


############
# Table A7 #
############

screenreg(
          list(
                mod.val.exclude,
                mod.glob.exclude,
                val.imp.exclude,
                glob.imp.exclude
               ),
          stars              = c(
                                 0.01,
                                 0.05,
                                 0.1
                                 ) ,
          custom.coef.names  = c(
                                 "Intercept",
                                 "Total trade (logged)",
                                 "Colonial ties",
                                 "Capabilities",
                                 "Africa",
                                 "Asia",
                                 "Failed mediation",
                                 "Ongoing mediation",
                                 "Intermediate exports (logged)",
                                 "Time",
                                 "Time$^{2}$",
                                 "Time$^{3}$",
                                 "Share of intermediate exports",
                                 "Intermediate imports (logged)",
                                 "Share of intermediate imports"
                                 ),
          custom.model.names = c(
                                 "Logged Value",
                                 "Global Share",
                                 "Logged Value",
                                 "Global Share"
                                 )
          )

############################################
# Omit regime type, GDPpc, and total trade #
############################################

mod.val.nocont <- glm(
                      majpow_med ~
                                   colonial +
                                   cinc + 
                                   africa +
                                   asia +
                                   num.fail +
                                   num.ongo +
                                   ln.exp + 
                                   t +
                                   t2 +
                                   t3,
                      data   = data.set,
                      family = binomial(link = "logit"),
                      x      = TRUE
                      )

mod.glob.nocont <- glm(
                       majpow_med ~
                                    colonial +
                                    cinc + 
                                    africa +
                                    asia +
                                    num.fail +
                                    num.ongo +
                                    expshare.glob + 
                                    t +
                                    t2 +
                                    t3,
                       data   = data.set,
                       family = binomial(link = "logit"),
                       x      = TRUE
                       )

val.imp.nocont <- glm(
                      majpow_med ~
                                   colonial +
                                   cinc + 
                                   africa +
                                   asia +
                                   num.fail +
                                   num.ongo +
                                   ln.exp + 
                                   ln.imp +
                                   t +
                                   t2 +
                                   t3,
                      data   = data.set,
                      family = binomial(link = "logit"),
                      x      = TRUE
                      )

glob.imp.nocont <- glm(
                       majpow_med ~
                                    colonial +
                                    cinc + 
                                    africa +
                                    asia +
                                    num.fail +
                                    num.ongo +
                                    expshare.glob + 
                                    impshare.glob +
                                    t +
                                    t2 +
                                    t3,
                       data   = data.set,
                       family = binomial(link = "logit"),
                       x      = TRUE
                       )

############
# Table A8 #
############

screenreg(
          list(
                mod.val.nocont,
                mod.glob.nocont,
                val.imp.nocont,
                glob.imp.nocont
               ),
          stars              = c(
                                 0.01,
                                 0.05,
                                 0.1
                                 ) ,
          custom.coef.names  = c(
                                 "Intercept",
                                 "Colonial ties",
                                 "Capabilities",
                                 "Africa",
                                 "Asia",
                                 "Failed mediation",
                                 "Ongoing mediation",
                                 "Intermediate exports (logged)",
                                 "Time",
                                 "Time$^{2}$",
                                 "Time$^{3}$",
                                 "Share of intermediate exports",
                                 "Intermediate imports (logged)",
                                 "Share of intermediate imports"
                                 ),
          custom.model.names = c(
                                 "Logged Value",
                                 "Global Share",
                                 "Logged Value",
                                 "Global Share"
                                 )
          )

###################################
# Omit all country-level controls #
###################################

mod.val.base <- glm(
                    majpow_med ~
                                 num.fail +
                                 num.ongo +
                                 ln.exp + 
                                 t +
                                 t2 +
                                 t3,
                    data   = data.set,
                    family = binomial(link = "logit"),
                    x      = TRUE
                    )


mod.glob.base <- glm(
                     majpow_med ~
                                  num.fail +
                                  num.ongo +
                                  expshare.glob + 
                                  t +
                                  t2 +
                                  t3,
                     data   = data.set,
                     family = binomial(link = "logit"),
                     x      = TRUE
                     )


val.imp.base <- glm(
                    majpow_med ~
                                 num.fail +
                                 num.ongo +
                                 ln.exp +
                                 ln.imp +
                                 t +
                                 t2 +
                                 t3,
                    data   = data.set,
                    family = binomial(link = "logit"),
                    x      = TRUE
                    )


glob.imp.base <- glm(
                     majpow_med ~
                                  num.fail +
                                  num.ongo +
                                  expshare.glob + 
                                  impshare.glob +
                                  t +
                                  t2 +
                                  t3,
                     data   = data.set,
                     family = binomial(link = "logit"),
                     x      = TRUE
                     )



############
# Table A9 #
############

screenreg(
          list(
                mod.val.base,
                mod.glob.base,
                val.imp.base,
                glob.imp.base
               ),
          stars              = c(
                                 0.01,
                                 0.05,
                                 0.1
                                 ) ,
          custom.coef.names  = c(
                                 "Intercept",
                                 "Failed mediation",
                                 "Ongoing mediation",
                                 "Intermediate exports (logged)",
                                 "Time",
                                 "Time$^{2}$",
                                 "Time$^{3}$",
                                 "Share of intermediate exports",
                                 "Intermediate imports (logged)",
                                 "Share of intermediate imports"
                                 ),
          custom.model.names = c(
                                 "Logged Value",
                                 "Global Share",
                                 "Logged Value",
                                 "Global Share"
                                 )
          )

